 use "C:\"your directory"\Italy for plot.dta" 

 *Figure 1
twoway (area stringencybis date, col(bluishgray)) (mspline new_deaths date, col(red) lpattern(dash)) (mspline new_cases date, col(blue) yaxis(2)) , legend(cols(3) label(1 "Stringency") label(2 "Deaths") label(3 "Cases") region(lwidth(none))) tlabel(01jan2020(93)01jan2022) xlabel(, labsize(small)) ylabel(, labsize(small)) ylabel(, axis(2) labsize(small)) ytitle("New deaths (per million) / Stringency", size(small)) ytitle("New cases (per million)", axis(2) size(small)) xtitle("") graphregion(color(white)) bgcolor(white)


use "C:\"your directory"\Italy mobility daily.dta" 

xtset codeprovince date

*Figure 2
twoway (mspline retail date, col(orange) lpattern(longdash)) (mspline transit date, col(midgreen) lpattern(shortdash)) (mspline workplaces date, col(midblue)) (mspline residential date, col(red) lpattern(dash)), legend(cols(4) label(1 "Retail") label(2 "Transit") label(3 "Workplaces") label(4 "Residential") region(lwidth(none))) tlabel(22159 22189 22220 22250 22281 22312 22340 22371 22401 22432 22462, format("%tdMon-YY")) xlabel(, labsize(small)) ylabel(, labsize(small)) ytitle("") xtitle("") graphregion(color(white)) bgcolor(white) yline(0, lcolor(black))


*Table A.1
xtpcse retail i.policy
xtpcse transit i.policy
xtpcse workplaces i.policy
xtpcse residential i.policy


*Table A.2 (producing the provincial coefficients mapped in Figure 4) (see below for the actual Maps code)
*b3=Ancona
xtpcse retail i.policy b3.codeprovince
xtpcse transit i.policy b3.codeprovince
xtpcse workplaces i.policy b3.codeprovince
xtpcse residential i.policy b3.codeprovince



*Figure 3
*Combomarginsplot
foreach var in retail transit workplaces residential {
      quietly xtpcse `var' i.policy
      margins policy, saving (`var'_1)
  }

combomarginsplot retail_1 transit_1 workplaces_1 residential_1, labels( "Retail" "Transit" "Workplaces" "Residential") plot3opts(lpattern(solid) msymbol(smcircle) mcolor(midblue) lcolor(midblue)) ci3opts(lcolor(midblue)) plot4opts(lpattern(dash) msymbol(smdiamond_hollow) mcolor(red) lcolor(red)) ci4opts(lcolor(red)) plot2opts(lpattern(shortdash) msymbol(smsquare) mcolor(midgreen) lcolor(midgreen)) ci2opts(lcolor(midgreen)) plot1opts(lpattern(longdash) msymbol(triangle) mcolor(orange) lcolor(orange)) ci1opts(lcolor(orange)) legend(cols(4) region(lwidth(none))) xlabel(1 "White" 2 "Yellow" 3 "Orange" 4 "Red", labsize(small)) ylabel(, labsize(small)) ytitle("") xtitle("") graphregion(color(white)) bgcolor(white) yline(0, lcolor(black)) title("")


*Table 1

xtpcse retail i.policy positivity fullypct service dipendenti gdppc1000 emprate soccapital EQI2021 logpopprov densityprov100 old i.dow, pairwise
xtpcse transit i.policy positivity fullypct service dipendenti gdppc1000 emprate soccapital EQI2021 logpopprov densityprov100 old i.dow, pairwise
xtpcse workplace i.policy positivity fullypct service dipendenti gdppc1000 emprate soccapital EQI2021 logpopprov densityprov100 old i.dow, pairwise
xtpcse residential i.policy positivity fullypct service dipendenti gdppc1000 emprate soccapital EQI2021 logpopprov densityprov100 old i.dow, pairwise


*Figure A.1

* A) First normalize dividing for 2sd (Gelman 2008)

foreach var in positivity fullypct service dipendenti gdppc1000 emprate soccapital EQI2021 logpopprov densityprov100 old {
      sum `var'
	  generate z`var' = `var'/(2*r(sd))
  }



* B) Run series of regressions and store estimates

foreach var in zpositivity zfullypct zservice zdipendenti zgdppc1000 zemprate zsoccapital zEQI2021 zlogpopprov zdensityprov100 zold {
      quietly xtpcse retail `var' i.policy i.dow
      estimates store retail`var'
  }


foreach var in zpositivity zfullypct zservice zdipendenti zgdppc1000 zemprate zsoccapital zEQI2021 zlogpopprov zdensityprov100 zold  {
      quietly xtpcse transit `var' i.policy i.dow
      estimates store transit`var'
  }

  
foreach var in zpositivity zfullypct zservice zdipendenti zgdppc1000 zemprate zsoccapital zEQI2021 zlogpopprov zdensityprov100 zold {
      quietly xtpcse workplace `var' i.policy i.dow
      estimates store workplace`var'
  }

  
foreach var in zpositivity zfullypct zservice zdipendenti zgdppc1000 zemprate zsoccapital zEQI2021 zlogpopprov zdensityprov100 zold  {
      quietly xtpcse residential `var' i.policy i.dow
      estimates store residential`var'
  }
  
  
* C) Plot coefficients

coefplot (retailzpositivity) (retailzfullypct) (retailzservice) (retailzdipendenti) (retailzgdppc1000) (retailzemprate) (retailzsoccapital) (retailzEQI2021) (retailzlogpopprov) (retailzdensityprov100) (retailzold), drop (_cons *.policy *.dow) graphregion(color(white)) bgcolor(white) legend(off) xtitle ("Retail") xline(0, lcolor(red)) ylabel(1 "Positivity" 2 "Vaccines" 3 "Service" 4 "Employees" 5 "GDP pc" 6 "Employment" 7 "Social capital" 8 "Gov quality" 9 "Population" 10 "Density" 11 "Old population", labsize(medsmall)) yla(, tlength(0)) 

coefplot (transitzpositivity) (transitzfullypct) (transitzservice) (transitzdipendenti) (transitzgdppc1000) (transitzemprate) (transitzsoccapital) (transitzEQI2021) (transitzlogpopprov) (transitzdensityprov100) (transitzold), drop (_cons *.policy *.dow) graphregion(color(white)) bgcolor(white) legend(off) xtitle ("Transit") xline(0, lcolor(red)) ylabel(1 "Positivity" 2 "Vaccines" 3 "Service" 4 "Employees" 5 "GDP pc" 6 "Employment" 7 "Social capital" 8 "Gov quality" 9 "Population" 10 "Density" 11 "Old population", labsize(medsmall)) yla(, tlength(0)) 

coefplot (workplacezpositivity) (workplacezfullypct) (workplacezservice) (workplacezdipendenti) (workplacezgdppc1000) (workplacezemprate) (workplacezsoccapital) (workplacezEQI2021) (workplacezlogpopprov) (workplacezdensityprov100) (workplacezold), drop (_cons *.policy *.dow) graphregion(color(white)) bgcolor(white) legend(off) xtitle ("Workplace") xline(0, lcolor(red)) ylabel(1 "Positivity" 2 "Vaccines" 3 "Service" 4 "Employees" 5 "GDP pc" 6 "Employment" 7 "Social capital" 8 "Gov quality" 9 "Population" 10 "Density" 11 "Old population", labsize(medsmall)) yla(, tlength(0)) 

coefplot (residentialzpositivity) (residentialzfullypct) (residentialzservice) (residentialzdipendenti) (residentialzgdppc1000) (residentialzemprate) (residentialzsoccapital) (residentialzEQI2021) (residentialzlogpopprov) (residentialzdensityprov100) (residentialzold), drop (_cons *.policy *.dow) graphregion(color(white)) bgcolor(white) legend(off) xtitle ("Residential") xline(0, lcolor(red)) ylabel(1 "Positivity" 2 "Vaccines" 3 "Service" 4 "Employees" 5 "GDP pc" 6 "Employment" 7 "Social capital" 8 "Gov quality" 9 "Population" 10 "Density" 11 "Old population", labsize(medsmall)) yla(, tlength(0)) 




*Table 2

xtpcse newcases l7.newcases l7.retail newtest l7.fullypct positivity i.dow, pairwise
xtpcse newcases l7.newcases l7.transit newtest l7.fullypct positivity i.dow, pairwise
xtpcse newcases l7.newcases l7.workplaces newtest l7.fullypct positivity i.dow, pairwise
xtpcse newcases l7.newcases l7.residential newtest l7.fullypct positivity i.dow, pairwise


*Table 3

xtpcse newcases l7.newcases l7.retail l7.i.policy newtest l7.fullypct positivity i.dow, pairwise
xtpcse newcases l7.newcases l7.transit l7.i.policy newtest l7.fullypct positivity i.dow, pairwise
xtpcse newcases l7.newcases l7.workplaces l7.i.policy newtest l7.fullypct positivity i.dow, pairwise
xtpcse newcases l7.newcases l7.residential l7.i.policy newtest l7.fullypct positivity i.dow, pairwise




*Figure 4
use "C:\Users\marco\Dropbox\a Covid\Mobility\ProvCM01012018_WGS84.dta"

grmap retail, fc(Blues2) ndsize(none) osize(none none none none none none none none none) cln(9) title(Retail and recreation, size(medium)) leg(off)
graph save "Graph" "A1 retail blu.gph", replace

grmap transit, fc(Blues2) ndsize(none) osize(none none none none none none none none none) cln(9) title(Transit stations, size(medium)) leg(off)
graph save "Graph" "A2 transit blu.gph", replace

grmap workplaces, fc(Blues2) ndsize(none) osize(none none none none none none none none none) cln(9) title(Workplaces, size(medium)) leg(off)
graph save "Graph" "A3 workplaces blu.gph", replace

ssc install palettes, replace
ssc install colrspace, replace
colorpalette Blues, select(2 2 3 4 5 6 7 8 9) nograph reverse
local colors `r(p)'
grmap residential, fc("`colors'") ndsize(none) osize(none none none none none none none none none) cln(9) title(Residential, size(medium)) leg(off)
graph save "Graph" "A4 residential blu.gph", replace

graph combine "A1 retail blu.gph" "A2 transit blu.gph" "A3 workplaces blu.gph" "A4 residential blu.gph", graphregion(color(white))
*Then adjust size 8" x 11"
